Uvod¶
Gibanje masne točke v polju sil v eni dimenziji opišemo z diferencialno enačbo drugega reda, z Newtonovim zakonom
\begin{equation*} m\, \frac{d^2x}{dt^2} = F \>. \end{equation*}
Enačba je seveda enakovredna sistemu enačb prvega reda $$ m \, \frac{dx}{dt} = p \;, \qquad \frac{dp}{dt} = F $$ in tako jo tudi rešujemo: kot sistem dveh enačb prvega reda.
Seveda morajo biti na voljo tudi ustrezni začetni pogoji, tipično $x(t=0)=x_0$ in $dx/dt=v(t=0)=v_0$. Splošnejše gre tu za sistem diferencialnih enačb drugega reda: $$ \frac{d^ny}{dx^n} = f(x,y,y',y'',...), $$ ki ga lahko prevedemo na sistem enačb prvega reda z uvedbo novih spremenljivk v slogu gibalne količine pri Netwonovi enačbi ($y'=v,y''=z,...$).
Z nekaj truda se da eksplicitno dokazati, mi pa lahko privzamemo, da so metode za reševanje enačb hoda (Runge-Kutta 4. reda, prediktor-korektor... ) neposredno uporabne za reševanje takšnih sistemov enačb in torej aplikabilne v poljubno dimenzijah, kar naj bi v principu zadovoljilo večino naših zahtev.
Obstaja še posebna kategorija tako imenovanih simplektičnih metod, za enačbe, kjer je $f$ le funkcija koordinat, $f(y)$, ki (približno) ohranjajo tudi Hamiltonian, torej energijo sistema. Najbolj znana metoda je Verlet/Stormer/Encke metoda, ki je globalno natančna do drugega reda in ki točno ohranja tudi vrtilno količino sistema (če je ta v danem problemu smiselna). Rešujemo torej za vsak diskretni korak $n$ velikosti $h$, $x_n=x_0+n \cdot h$:
$$ \frac{d^2y}{dx^2} = f(y) $$
in pri diskretizaciji dobimo recept za korak $y_n$ in $v_n=y'_n$:
$$ y_{n+1} = y_n + h \cdot v_n + \frac{h^2}{2} \cdot f(y_n) \\$$ $$ v_{n+1} = v_n + \frac{h}{2} \cdot \left[ f(y_n) + f(y_{n+1}) \right]. $$ Alternativno lahko to shemo zapišemo tudi s pomočjo dodatnih vmesnih točk in preskakujemo med lego in hitrostjo z zamikom h/2 (od tod angleško ime 'leapfrog' za ta zapis): $$ y_{n+1} = y_n + h \cdot v_{n+1/2} \ $$ $$ v_{n+3/2} = v_{n+1/2} + h \cdot f(y_{n+1}). $$ V še enem drugačnem zapisu je metoda poznana tudi kot metoda *Središčne razlike* (Central Difference Method, CDM), če nas hitrost ne zanima: $$ y_{n+1} - 2 y_n + y_{n-1} = h^2 \cdot f(y_n), $$
kjer prvo točko $y_1$ izračunamo po originalni shemi. Metodo CDM lahko uporabljamo tudi za primere, ko je f tudi funkcija 'časa' x, f(x,y), le da tu simplektičnost ni zagotovljena (in tudi verjetno ne relevantna). Za simplektične metode višjih redov je na voljo na primer Forest-Ruth metoda ali Position Extended Forest-Ruth Like (PEFRL) metoda, ki sta obe globalno četrtega reda in enostavni za implementacijo.
Naloga¶
Čim več metod uporabi za izračun nihanja matematičnega nihala z začetnim pogojem $\theta(0)= \theta_0 = 1$, $\dot{\theta}(0)=0$. Poišči korak, ki zadošča za natančnost na 3 mesta. Primerjaj tudi periodično stabilnost shem: pusti, naj teče račun čez 10 ali 20 nihajev in poglej, kako se amplitude nihajev sistematično kvarijo. Pomagaš si lahko tudi tako, da občasno izračunaš energijo $E \propto 1-\cos \theta + \frac{\dot{\theta}^2 }{2 \omega_0^2} $. Nariši tudi ustrezne fazne portrete!.
Dodatno: Razišči še resonančno krivuljo vzbujenega dušenega matematičnega nihala \begin{equation*} \frac{d^2x}{dt^2} + \beta\frac{d^2x}{dt^2}+ \sin x = v\cos\omega_0t\>, \end{equation*} kjer je $\beta$ koeficient dušenja, $v$ in $\omega_0$ pa amplituda in frekvenca vzbujanja. Opazuj obnašanje odklonov in hitrosti nihala pri dušenju $\beta=0.5$, vzbujevalni frekvenci $\omega_0=2/3$ in amplitudo vzbujanja na območju $0.5 < v < 1.5$. Poskusi opaziti histerezno obnašanje resonančne krivulje pri velikih amplitudah vzbujanja
Matematično nihalo¶
Nihanje matematičnega nihala lahko opišemo z diferencialno enačbo:
$$ \frac{d^2\theta}{dt^2} + \omega_0 \sin(\theta) = 0. $$
Če so odmiki majhni, lahko enačbo linearno aproksimiramo in dobimo harmonične rešitve kot:
$$ \theta (t) = \theta_0 \sin(\omega_0 t + \phi_0), $$
kjer sta $\theta_0$ in $\phi_0$ odvisna od začetnih pogojev. V splošnem pa zgodba ni tako preprosta in se moramo za nihajni čas in analitično rešitev časovnega poteka zateči k specialnim funkcijam. Nihajni čas lahko izračunamo s pomočjo eliptičnega integrala prve vrste, in sicer:
$$ T_0 = \frac{4}{\omega_0} \int_0^{\pi / 2} \frac{d \varphi}{\sqrt{1 - \sin^2(\theta_0 /2) \sin^2 \varphi}} = \frac{4}{\omega_0} K(\sin^2(\theta_0 / 2)), $$
kjer $\theta_0$ predstavlja amplitudo nihanja. Časovni potek odmika pa je podan s pomočjo nepopolne eliptične funkcije prve vrste oz. Jacobijeve amplitudne funkcije $sn(u;m)$. Za amplitudo nihanja sledi
$$ \theta(t) = 2 \arcsin \left( \sin \frac{\theta_0}{2} \cdot sn(\omega_0 (T_0 / 4 - t), \sin^2(\theta_0 / 2)) \right).$$
Ilustracija rešitev DE¶
Za začetek narišimo nekaj referenčnih grafov s pomočjo analitične rešitve.
Nadaljujemo s fazno sliko rešitve DE. V teoriji bi enake rezultate dobili z odvajanjem analitične rešitve. A se tokrat kot predpripava za nadaljne raziskovanje odločimo rešitve pridobiti numerično s simplektično metodo Verlet.
Opazujemo parametrično krivuljo v faznem prostoru. Za majhne začetne hitrosti so rešitve v faznem prostoru krožnice, za večje začetne hitrosti pa se popačijo v obliko podobno očesu.
Za zabavo lahko skiciramo še analitičen graf odvisnosti časa $T_0$ od začetne lege. S pomočjo funkcije scipy.signal.find_peaks poiščemo in izmerimo intervale $T_0$. Ugotovimo, da se točke v okviru napake (ki je v tem primeru kar enaka dvema dolžinama integracijskega koraka = 0.2) ujemajo z analitično napovedjo. Ugotovimo tudi, da časa $T_0$ numerično ne bomo znali določiti zelo natančno, če ne bo naš korak res zelo majhen.
Primerjava integracijskih metod¶
Sedaj primerjamo med sabo različne numerične metode za integracijo problema. Od prej poznamo že Eulerjevo metodo, Heunovo metodo, Runge-Kutta 4. reda, le-tem pa dodamo še 2 simplektični; Verletovo metodo in PEFRL.
S funkcijo find_peaks sedaj preverimo še kaj se dogaja z amplitudo skozi dolg čas. Ugotovimo, da se napaka sistematično veča. Če pogledamo sedaj še kaj se dogaja z energijo sistema skozi čas, ugotovimo, da se le ta sistematično spreminja. Le pri simplektičnih metodah sicer niha, ampak skozi daljši čas ostaja konstantna - kar je tudi naš cilj.
Poglejmo si še, kako se simplektične in ostale metode razlikujejo na faznem portretu. Kot vidimo na spodnjem grafu, se rešitev s pomočjo metode Runge-Kutta po dolgem času vrača v ravnovesno lego. Medtem ko rešitev s simplektično metodo Verlet ohranja "orbito" oz. energijo.
Dodatek: Resonančna krivulja vzbujenega dušenega matematičnega nihala¶
V dodatku analiziramo vzbujano dušeno matematično nihalo. Nihalo opišemo s tremi parametri; dušilnim faktorjem $\beta$, vzbujalo amplitudo $v$ in vzbujalno frevenco $\omega$. Začetni pogoji so enaki kot pri prejšnji nalogi; $\theta = 1$ in $\dot \theta = 0$. Opazimo, da se sistem za nekatere vrednosti parametrov obnaša kaotično. Pri zmernih vrednosti $v$, se sistem obnaša mirno, pričakovano. Za večje vrednosti pa pobezlja. Podobno lahko opazimo tudi pri drugačnem spreminjanju parametrov. Na koncu izrišemo še resonančno krivuljo, ki je stabilna le nekje do $v \approx 0.8 $.
Zaključek¶
Tekom te naloge smo primerjali različne metode integracije za reševanje diferencialnih enačb matematičnega nihala. Ugotovili smo, da simplektične metode, kot sta Verletova in PEFRL, najboljše ohranjajo energijo sistema skozi čas, medtem ko klasične metode, kot je Runge-Kutta, lahko vodijo do sistematičnega kopičenja napak. Analizirali smo tudi resonančno krivuljo vzbujenega dušenega nihala in opazili kaotično obnašanje pri večjih amplitudah vzbujanja. Izkušnje kažejo, da je izbira metode ključna za natančnost in stabilnost dolgoročnih simulacij.
Luka Skeledžija, Github source 🔗, 2024